Stiff ordinary differential equations

Consider the Lotka-Voterra system from HW1: \begin{align*} \begin{cases} u_1'(t) = u_1(t) - u_1(t) u_2(t),\\ u_2'(t) = u_1(t) u_2(t) - u_2(t). \end{cases} \end{align*} For $u_1(0) = 5, u_2(0) = 0.8$.

We add a third variable to this model that gives some interesting dynamics:

\begin{align*} \begin{cases} u_1'(t) = u_1(t) - u_1(t) u_2(t) + \kappa_1 u_1(t)u_2(t)u_3(t),\\ u_2'(t) = u_1(t) u_2(t) - u_2(t)+ \kappa_2 u_1(t)u_2(t)u_3(t),\\ u_3'(t) = -\kappa_3 u_1(t)u_2(t)u_3(t). \end{cases} \end{align*}

We suppose that $\kappa_1 + \kappa_2 = \kappa_3$.

This forces the third variable $u_3(t)$ to be severely dampled with $u_1(t)u_2(t)$ is large.

Choose $\kappa_1 = \varkappa/2 = \kappa_2$, $\kappa_3 = \varkappa$.

We might think that because $u_3(t)$ is damped, the approximation of it is irrelevant in the long term. But if we start with a different initial condition for $u_3(0)$, we get a different limit cycle:

We see that the amplitude of the oscillations has increased. So, if we accumulate error in our approximation of $u_3(t)$, it will necessarily effect the quality of the solution for all time. This is different than for the equation

\begin{align*} u'(t) = \lambda(u(t) - h(t)) + h'(t), \end{align*}

where initial errors will be damped out overtime and the limiting solution $u(t) = h(t)$ does not depend on the choice of initial condition.

Forward Euler fails on this problem with this time step due to instability. But if we roll back the $\varkappa$ we can start to see what is happening.

Compare this to the solution where $u_3(t) = 0$:

So, capturing the initial transient is extremely important in determining the behavior of the solution. But we have to shrink the time step just to capture this transient. This is stiffness --- severe differences in timescales.